{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "H9h3CaW8kj1y"
   },
   "source": [
    "**Problem setup**  \n",
    "  \n",
    "We are going to solve the non-linear Schrödinger equation given by  \n",
    "$i h_t + \\frac{1}{2} h_{xx} + |h|^2h = 0$  \n",
    "  \n",
    "with periodic boundary conditions as  \n",
    "$x \\in [-5,5], \\quad t \\in [0, \\pi/2]$  \n",
    "$h(t, -5) = h(t,5)$  \n",
    "$h_x(t, -5) = h_x(t,5)$  \n",
    "  \n",
    "and initial condition equal to  \n",
    "$h(0,x) = 2 sech(x)$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "f7Fo1dqAmMt-"
   },
   "source": [
    "Deepxde only uses real numbers, so we need to explicitly split the real and imaginary parts of the complex PDE.  \n",
    "  \n",
    "In place of the single residual  \n",
    "$f = ih_t + \\frac{1}{2} h_{xx} +|h|^2 h$  \n",
    "  \n",
    "we get the two (real valued) residuals  \n",
    "$f_{\\mathcal{R}} = u_t + \\frac{1}{2} v_{xx} + (u^2 + v^2)v$  \n",
    "$f_{\\mathcal{I}} = v_t - \\frac{1}{2} u_{xx} - (u^2 + v^2)u$  \n",
    "  \n",
    "where u(x,t) and v(x,t) denote respectively the real and the imaginary part of h.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'2.3.0'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[PhysicalDevice(name='/physical_device:GPU:0', device_type='GPU')]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf.config.list_physical_devices('GPU')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import torch\n",
    "\n",
    "torch.cuda.is_available()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 6652,
     "status": "ok",
     "timestamp": 1640283904989,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "YjOFqj77nl_R",
    "outputId": "aa0997be-fdbe-4a21-8def-580cf5e20c6c"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using backend: tensorflow.compat.v1\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\tensorflow\\python\\compat\\v2_compat.py:96: disable_resource_variables (from tensorflow.python.ops.variable_scope) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "non-resource variables are not supported in the long term\n",
      "WARNING:tensorflow:From G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\deepxde\\nn\\initializers.py:118: The name tf.keras.initializers.he_normal is deprecated. Please use tf.compat.v1.keras.initializers.he_normal instead.\n",
      "\n",
      "WARNING:tensorflow:From G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\deepxde\\nn\\initializers.py:119: The name tf.keras.initializers.he_uniform is deprecated. Please use tf.compat.v1.keras.initializers.he_uniform instead.\n",
      "\n",
      "WARNING:tensorflow:From G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\deepxde\\nn\\initializers.py:120: The name tf.keras.initializers.lecun_normal is deprecated. Please use tf.compat.v1.keras.initializers.lecun_normal instead.\n",
      "\n",
      "WARNING:tensorflow:From G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\deepxde\\nn\\initializers.py:121: The name tf.keras.initializers.lecun_uniform is deprecated. Please use tf.compat.v1.keras.initializers.lecun_uniform instead.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "import deepxde as dde\n",
    "\n",
    "# For plotting\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.interpolate import griddata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "executionInfo": {
     "elapsed": 337,
     "status": "ok",
     "timestamp": 1640283910535,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "Oj-6F7RuobUn"
   },
   "outputs": [],
   "source": [
    "x_lower = -5\n",
    "x_upper = 5\n",
    "t_lower = 0\n",
    "t_upper = np.pi / 2\n",
    "\n",
    "# Creation of the 2D domain (for plotting and input)\n",
    "x = np.linspace(x_lower, x_upper, 256)\n",
    "t = np.linspace(t_lower, t_upper, 201)\n",
    "X, T = np.meshgrid(x, t)\n",
    "\n",
    "# The whole domain flattened\n",
    "X_star = np.hstack((X.flatten()[:, None], T.flatten()[:, None]))\n",
    "\n",
    "# Space and time domains/geometry (for the deepxde model)\n",
    "space_domain = dde.geometry.Interval(x_lower, x_upper)\n",
    "time_domain = dde.geometry.TimeDomain(t_lower, t_upper)\n",
    "geomtime = dde.geometry.GeometryXTime(space_domain, time_domain)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "executionInfo": {
     "elapsed": 283,
     "status": "ok",
     "timestamp": 1640283914251,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "dbx1YulhoORs"
   },
   "outputs": [],
   "source": [
    "# The \"physics-informed\" part of the loss\n",
    "\n",
    "\n",
    "def pde(x, y):\n",
    "    \"\"\"\n",
    "    INPUTS:\n",
    "        x: x[:,0] is x-coordinate\n",
    "           x[:,1] is t-coordinate\n",
    "        y: Network output, in this case:\n",
    "            y[:,0] is u(x,t) the real part\n",
    "            y[:,1] is v(x,t) the imaginary part\n",
    "    OUTPUT:\n",
    "        The pde in standard form i.e. something that must be zero\n",
    "    \"\"\"\n",
    "\n",
    "    u = y[:, 0:1]\n",
    "    v = y[:, 1:2]\n",
    "\n",
    "    # In 'jacobian', i is the output component and j is the input component\n",
    "    u_t = dde.grad.jacobian(y, x, i=0, j=1)\n",
    "    v_t = dde.grad.jacobian(y, x, i=1, j=1)\n",
    "\n",
    "    u_x = dde.grad.jacobian(y, x, i=0, j=0)\n",
    "    v_x = dde.grad.jacobian(y, x, i=1, j=0)\n",
    "\n",
    "    # In 'hessian', i and j are both input components. (The Hessian could be in principle something like d^2y/dxdt, d^2y/d^2x etc)\n",
    "    # The output component is selected by \"component\"\n",
    "    u_xx = dde.grad.hessian(y, x, component=0, i=0, j=0)\n",
    "    v_xx = dde.grad.hessian(y, x, component=1, i=0, j=0)\n",
    "\n",
    "    f_u = u_t + 0.5 * v_xx + (u ** 2 + v ** 2) * v\n",
    "    f_v = v_t - 0.5 * u_xx - (u ** 2 + v ** 2) * u\n",
    "\n",
    "    return [f_u, f_v]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "executionInfo": {
     "elapsed": 323,
     "status": "ok",
     "timestamp": 1640283979042,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "OisdDqf0oSLx"
   },
   "outputs": [],
   "source": [
    "# Boundary and Initial conditions\n",
    "\n",
    "# Periodic Boundary conditions\n",
    "bc_u_0 = dde.PeriodicBC(\n",
    "    geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=0, component=0\n",
    ")\n",
    "bc_u_1 = dde.PeriodicBC(\n",
    "    geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=1, component=0\n",
    ")\n",
    "bc_v_0 = dde.PeriodicBC(\n",
    "    geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=0, component=1\n",
    ")\n",
    "bc_v_1 = dde.PeriodicBC(\n",
    "    geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=1, component=1\n",
    ")\n",
    "\n",
    "# Initial conditions\n",
    "def init_cond_u(x):\n",
    "    \"2 sech(x)\"\n",
    "    return 2 / np.cosh(x[:, 0:1])\n",
    "\n",
    "\n",
    "def init_cond_v(x):\n",
    "    return 0\n",
    "\n",
    "\n",
    "ic_u = dde.IC(geomtime, init_cond_u, lambda _, on_initial: on_initial, component=0)\n",
    "ic_v = dde.IC(geomtime, init_cond_v, lambda _, on_initial: on_initial, component=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "executionInfo": {
     "elapsed": 259,
     "status": "ok",
     "timestamp": 1640283983373,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "u-f8sfj2oWgy"
   },
   "outputs": [],
   "source": [
    "data = dde.data.TimePDE(\n",
    "    geomtime,\n",
    "    pde,\n",
    "    [bc_u_0, bc_u_1, bc_v_0, bc_v_1, ic_u, ic_v],\n",
    "    num_domain=10000,\n",
    "    num_boundary=20,\n",
    "    num_initial=200,\n",
    "    train_distribution=\"pseudo\",\n",
    ")\n",
    "\n",
    "# Network architecture\n",
    "net = dde.nn.FNN([2] + [100] * 4 + [2], \"tanh\", \"Glorot normal\")\n",
    "\n",
    "model = dde.Model(data, net)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Hu5Tyq32DezT"
   },
   "source": [
    "Adam optimization.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 439512,
     "status": "ok",
     "timestamp": 1640284427354,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "BR_s0D-_oZYz",
    "outputId": "9071ee58-c123-427b-fe32-cc8fa7c721f4",
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compiling model...\n",
      "Building feed-forward neural network...\n",
      "WARNING:tensorflow:From G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\deepxde\\nn\\tensorflow_compat_v1\\fnn.py:103: dense (from tensorflow.python.keras.legacy_tf_layers.core) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Use keras.layers.Dense instead.\n",
      "WARNING:tensorflow:From G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\tensorflow\\python\\keras\\legacy_tf_layers\\core.py:187: Layer.apply (from tensorflow.python.keras.engine.base_layer_v1) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Please use `layer.__call__` method instead.\n",
      "'build' took 0.058478 s\n",
      "\n",
      "'compile' took 1.530092 s\n",
      "\n",
      "Initializing variables...\n",
      "Training model...\n",
      "\n",
      "Step      Train loss                                                                          Test loss                                                                           Test metric\n",
      "0         [5.93e-03, 5.88e-04, 3.36e-01, 1.86e-06, 2.36e-02, 3.18e-05, 8.05e-01, 5.97e-03]    [5.93e-03, 5.88e-04, 3.36e-01, 1.86e-06, 2.36e-02, 3.18e-05, 8.05e-01, 5.97e-03]    []  \n",
      "1000      [1.02e-02, 1.13e-02, 1.49e-05, 3.22e-04, 3.67e-06, 1.40e-06, 1.99e-02, 2.25e-03]    [1.02e-02, 1.13e-02, 1.49e-05, 3.22e-04, 3.67e-06, 1.40e-06, 1.99e-02, 2.25e-03]    []  \n",
      "2000      [8.70e-03, 6.13e-03, 3.60e-04, 5.01e-05, 3.94e-04, 3.31e-05, 1.48e-02, 7.26e-04]    [8.70e-03, 6.13e-03, 3.60e-04, 5.01e-05, 3.94e-04, 3.31e-05, 1.48e-02, 7.26e-04]    []  \n",
      "3000      [6.84e-03, 3.50e-03, 3.22e-04, 1.54e-05, 2.44e-04, 2.16e-05, 1.03e-02, 2.92e-04]    [6.84e-03, 3.50e-03, 3.22e-04, 1.54e-05, 2.44e-04, 2.16e-05, 1.03e-02, 2.92e-04]    []  \n",
      "4000      [5.47e-03, 5.06e-03, 3.34e-05, 1.01e-05, 3.76e-05, 2.52e-05, 5.88e-03, 1.68e-04]    [5.47e-03, 5.06e-03, 3.34e-05, 1.01e-05, 3.76e-05, 2.52e-05, 5.88e-03, 1.68e-04]    []  \n",
      "5000      [2.17e-03, 2.32e-03, 4.12e-06, 8.06e-06, 3.13e-06, 2.49e-05, 4.02e-03, 2.43e-05]    [2.17e-03, 2.32e-03, 4.12e-06, 8.06e-06, 3.13e-06, 2.49e-05, 4.02e-03, 2.43e-05]    []  \n",
      "6000      [1.72e-03, 1.75e-03, 1.70e-06, 8.99e-06, 1.13e-06, 1.18e-05, 2.75e-03, 2.73e-05]    [1.72e-03, 1.75e-03, 1.70e-06, 8.99e-06, 1.13e-06, 1.18e-05, 2.75e-03, 2.73e-05]    []  \n",
      "7000      [1.37e-03, 1.34e-03, 1.45e-06, 9.40e-06, 1.05e-06, 7.21e-06, 1.98e-03, 2.09e-05]    [1.37e-03, 1.34e-03, 1.45e-06, 9.40e-06, 1.05e-06, 7.21e-06, 1.98e-03, 2.09e-05]    []  \n",
      "8000      [8.91e-04, 1.11e-03, 2.12e-06, 9.13e-06, 1.70e-06, 5.34e-06, 1.40e-03, 1.42e-05]    [8.91e-04, 1.11e-03, 2.12e-06, 9.13e-06, 1.70e-06, 5.34e-06, 1.40e-03, 1.42e-05]    []  \n",
      "9000      [1.41e-03, 6.76e-04, 1.26e-04, 6.35e-06, 7.27e-06, 6.37e-06, 1.13e-03, 1.22e-05]    [1.41e-03, 6.76e-04, 1.26e-04, 6.35e-06, 7.27e-06, 6.37e-06, 1.13e-03, 1.22e-05]    []  \n",
      "10000     [6.76e-04, 5.85e-04, 3.17e-06, 4.71e-06, 2.34e-06, 4.80e-06, 7.04e-04, 8.69e-06]    [6.76e-04, 5.85e-04, 3.17e-06, 4.71e-06, 2.34e-06, 4.80e-06, 7.04e-04, 8.69e-06]    []  \n",
      "\n",
      "Best model at step 10000:\n",
      "  train loss: 1.99e-03\n",
      "  test loss: 1.99e-03\n",
      "  test metric: []\n",
      "\n",
      "'train' took 218.871333 s\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(<deepxde.model.LossHistory at 0x22b287067c0>,\n",
       " <deepxde.model.TrainState at 0x22b7cf2ea00>)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# To employ a GPU accelerated system is highly encouraged.\n",
    "\n",
    "model.compile(\"adam\", lr=1e-3, loss=\"MSE\")\n",
    "model.train(epochs=10000, display_every=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xbAdby9zDosK"
   },
   "source": [
    "L-BFGS optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 438788,
     "status": "ok",
     "timestamp": 1640284883243,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "SVSKBTp6qRX8",
    "outputId": "f7ec1668-f24a-414e-c0bd-c207bb208505",
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compiling model...\n",
      "'compile' took 0.389958 s\n",
      "\n",
      "Training model...\n",
      "\n",
      "Step      Train loss                                                                          Test loss                                                                           Test metric\n",
      "10000     [1.34e-02, 6.86e-03, 1.69e-03, 1.87e-06, 8.76e-04, 8.16e-06, 1.24e-03, 7.40e-04]    [1.34e-02, 6.86e-03, 1.69e-03, 1.87e-06, 8.76e-04, 8.16e-06, 1.24e-03, 7.40e-04]    []  \n",
      "11000     [2.11e-05, 2.97e-05, 1.69e-07, 2.83e-07, 5.19e-08, 2.62e-07, 6.46e-06, 2.25e-07]                                                                                            \n",
      "12000     [6.28e-06, 7.91e-06, 3.62e-08, 3.12e-07, 2.07e-08, 2.97e-07, 3.76e-06, 1.89e-07]                                                                                            \n",
      "13000     [2.88e-06, 3.82e-06, 3.54e-08, 3.47e-07, 1.33e-08, 6.02e-08, 3.02e-06, 9.07e-08]                                                                                            \n",
      "14000     [1.87e-06, 2.25e-06, 4.49e-08, 1.91e-07, 1.49e-08, 1.49e-08, 2.66e-06, 7.09e-08]                                                                                            \n",
      "15000     [1.35e-06, 1.57e-06, 3.13e-08, 9.89e-08, 7.99e-09, 1.15e-08, 2.24e-06, 7.74e-08]                                                                                            \n",
      "INFO:tensorflow:Optimization terminated with:\n",
      "  Message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH\n",
      "  Objective function value: 0.000005\n",
      "  Number of iterations: 4748\n",
      "  Number of functions evaluations: 5051\n",
      "15051     [1.33e-06, 1.56e-06, 3.24e-08, 9.76e-08, 6.84e-09, 1.34e-08, 2.23e-06, 8.24e-08]    [1.33e-06, 1.56e-06, 3.24e-08, 9.76e-08, 6.84e-09, 1.34e-08, 2.23e-06, 8.24e-08]    []  \n",
      "\n",
      "Best model at step 15051:\n",
      "  train loss: 5.36e-06\n",
      "  test loss: 5.36e-06\n",
      "  test metric: []\n",
      "\n",
      "'train' took 173.443512 s\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(<deepxde.model.LossHistory at 0x1d73afd77c0>,\n",
       " <deepxde.model.TrainState at 0x1d73afd7430>)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dde.optimizers.config.set_LBFGS_options(\n",
    "    maxcor=50,\n",
    "    ftol=1.0 * np.finfo(float).eps,\n",
    "    gtol=1e-08,\n",
    "    maxiter=10000,\n",
    "    maxfun=10000,\n",
    "    maxls=50,\n",
    ")\n",
    "model.compile(\"L-BFGS\")\n",
    "model.train()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "u9EXdIk0FjYj"
   },
   "source": [
    "Final results.  \n",
    "The reference solution and further information can be found in [this paper](https://arxiv.org/abs/1711.10561) from Raissi, Karniadakis, Perdikaris.  \n",
    "The test data can be got [here](https://github.com/maziarraissi/PINNs/blob/master/main/Data/NLS.mat)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 281
    },
    "executionInfo": {
     "elapsed": 6595,
     "status": "ok",
     "timestamp": 1640284904382,
     "user": {
      "displayName": "Federico Magnani",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "09090763345999242901"
     },
     "user_tz": -60
    },
    "id": "aYl4TD96FrLV",
    "outputId": "30c43f0d-1147-40a7-f0c9-49d6e6466f65"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEICAYAAABWJCMKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABTf0lEQVR4nO29ebQ9SVXn+9mZ5/7qx1gFMkqBIAsbFQEFEZCnFIpMLeWAE6gt2hR0P9C2sdXWtYQWV7fa0uLUIiDa2tg4AQqWIDwU9GHJtKCgAH1lAVIUWpagFBT1+92Tud8fGREZmRmZGWfMPPfGd627Tt4Yd8bw3RE7hhRVJSEhISHh9CGbWoCEhISEhGmQFEBCQkLCKUVSAAkJCQmnFEkBJCQkJJxSJAWQkJCQcEqRFEBCQkLCKUVSAAkJe4CI/JmI/Nup5UhI8JEUQMKphIh8SEQ+IyKfEpG/F5FfF5Fb7ynv7xKRv9hHXgkJQ0gKIOE042tV9dbAA4EvBv7ztOIkJOwXSQEknHqo6t8Dr6NSBIjIQ0XkLSLyzyLybhF5pA1rRu/XiMiNIvJBEXmKcX+uiPxvL9w9RURFZOHnJSKfD7wQeJiZffyzcX+8iLzPpPtREfmB3b51QkJSAAkJiMjFwOOAq0XkbsAfAT8B3B74AeD3ReSOInIr4OeBx6nqbYCHA+9aJS9VfT/wDOAvVfXWqnqR8fpV4Okm3fsBb9z4xRISRpAUQMJpxqtE5EbgI8D1wHOAbwcuV9XLVbVU1dcDbwceb+KUwP1E5Baq+jFVvWpLshwDXyAit1XVT6jqO7eUbkJCL5ICSDjN+Doz4n4kcF/gDsDnAN9kzD//bEw0jwDuqqqfBr6FagT/MRH5IxG575Zk+UYqJfNhEXmTiDxsS+kmJPQiKYCEUw9VfRPw68DPUM0GflNVL/L+bqWqP2nCvk5VHw3cFfgA8GKTzKeBW3rJ3mUoy4AMb1PVS4E7Aa8Cfmezt0pIGEdSAAkJFV4APBr4C+BrReQxIpKLyFkReaSIXCwidxaRJ5q1gHPAp4DCxH8X8BUicg8RuZDhHUX/AFwsImcAROSMiDxFRC5U1WPgk166CQk7Q1IACQmAqv4j8BvAfwAuBX4E+EeqGcF/ouorGfBs4Drg48BXAv/exH898NvAlcA7gNcMZPdG4Crg70XkBuP2HcCHROSTVCamb9/e2yUkhCHpgzAJCQkJpxNpBpCQkJBwSrEYD7I7iMiHgBup7J1LVX3wlPIkJCQknCZMqgAMLlHVG8aDJSQkJCRsE8kElJCQkHBKMekisIh8EPgE1b7oX1HVFwXCXAZcBiBnjx509uI77FfIhISEhAPHZ67+2A2qese2+9QK4LNV9ToRuRPweuBZqvrmvvC3vM9n633+x/fsT8CEhISEE4Arn/gT7witsU5qAlLV68zv9cArgYdMKU9CQkLCacJkCkBEbiUit7HPwNcA751KnoSEhITThil3Ad0ZeKWIWDl+S1VfOxxFEUkH104jVGVqERISThwmUwCqeg3wgKnyTzgsJMV/spAU+jwwh3MACQkJpwxJoYexb8V4cAogNZyEk4g0Ik6A/fPbQSkAAbId95My6ZeECZAGNtPjNCrhg1IA+8CuFUxCwhyRBj6nUwkfnAI4jZWUcFg4xJFkGvjsDnNWroelAASypAASdoRyS8SdBimHg30o6zkr11EFICK/qarfMea2L6TOlbAr5F7bOsRRfEITMQr9tPNJzAzgC/1/RCQHHrQbcYYhaKOTxuB0V2/C2tgRMWxrlpEwjhBXJMXeRK8CEJH/TPVd1FuY75RCtRHnPNC5tXNfSCageSIRWxxWHcDsEvORZI/YQvmfpLbeqwBU9b+JyE8BL1HV796jTL0QgTwrpxZjEsx95HIaFPNJ6viHgjmW+a6U+BQ9aNAEpKqliMzmugZBWZwCBRBs9KeAYDfBPohiTqP3OWGXg5NDGVjMUVHFIGYN4AoR+VJVfdvOpYlALidfAeRrtqW5zxK2hZLuex4KURw60uCkhl8WUw8O1u37MQrgEuDpIvJh4NNU6wCqqvdfK8cNkIlyJi/2ne1WcIgjhLnKnA/4zVXmk4J1Byer4GAGMiNihgYqu8Q6bT9GATxudVF2AwHOZMupxVgJrlImatOlbvDJhw1k3nfjXxVJUSSsgnXay9BAZSf5rdGkRxWAqn4YwHy28ezqWWwPmShn890qgKmIYduEWb/HdkxmuyyXScp8bPS2ieKcGNkpMJNOgV23iU05YCczABF5IvB84LOB64HPAd5P63zAPiAot8iP95LXtklp08ot1pRnk0bbKIPI7LelyCZTxCbfTPZnapzr+sVc5ZoO67WJTdryOnFXiRNjAnoe8FDgDar6xSJyCfBtK0u1BeRScqvFubXibp3Q1yDWWHKMJfs4GYrod4+RbxVFtGoZ+XKuSz6bKqC5kF625U2Bm77XNmcVUy+Y7hvrDN5i+85Y3x7rDzEK4FhV/0lEMhHJVPVPzfmAvSMT5db5egrAxybKoIj4jPJopYz4D+UxFLev0YQawVCjHGp82UiD8vPKW6PosY5g70xZRXG0y2PbxLkPrEPOsYQcS7ax5TYm66qKYhNlMMe6XmcAMtYvYvpDiBdiZIlRAP8sIrcG/hx4mYhcD0yyEptLyYWLm1aKs+ootBgJP1SoQ3HDJNwNH6rIkEIIVrjnNhZnESBbK6NP3O3GmUm4TJ3pxOuU7Xf2O3uo0Y/V1TZnMrvEtogphkxDBNqXf4i8h/JYJe2+tPoURh6xNhWrGOcya+vDJgPCvsFguP8H+uUWZgCXAjcD/wF4CnAh8OMR8baOnJIL88+sFXezUXfXrwgU7Fi4oH/Lza+wIkDOQT9vp5F9j0y14293CTTf1euIAVksCQzPGLx3bMgvrXCRiixYtnFuY/JtC5uQU9utYfpqEGy3vFYdYft5teP6BO/nG0PkPoE38+i+b5vsg2FCeQaUzdg5oFWV7zrnisYGiW2MkXBwIDgyYIwZKMa0+5hdQJ8WkbsADwE+DrxOVf9pNOUdYCEln5V/KugXIuQQ+kaZIQUQHKW2wvmV4qfRLvymX2bcusoh6BZQBI2RtkgdTqsGHSL7MfNRV8mIkyFE8n75DIVruGnIrV8JjhH/0Dvtaj953w2SMWSfRZJuiFQzUaekbdwyRPD+a2u/bH46fv42bIjkh/x8MrXphtza7u3884FyyQMEH1JYYzOMbc0aYkg2PBuPHyy2ua3hNzJI3HgGICL/Fvgx4I1UTesXROTHVfWlY3G3jZyCi/JPR4fv09RtEoduITcLNlzg7Xg1ifvhxaVh4xbSJWXrVmrm0iy9zukq3f5KOTgSsTKUKnUeDTfz7JF5m8RLxMm/LHPn1ibxUHp+vvZXVToN0g+nQeVg/Fpx/PDt9CxiFUDpSDIqeEMB+ETSVgwNP8+tTT7iublf1KXn+/n+lVvZjSvKwpCiJcccJXM7nGz+pRPMtqRM1JGnn17bLZeyQ/IZWj/7+XrPvl/lVqfbJu0qj/5ZS1MJhZRAP8lv40aBvv4Xbx1oDyYlaoAZMvX6vOHnOTYwjjEB/Sfgi+2oX0Q+C3gLsH8FIMpFWbwJaB0zzdhIvXBE19XQvl9N9uLSsiTfUArWZGNHdJRktlJtA9Yc2g1Ws+CmtHYjOdbcPS81d/n7hA5wXPrh6jTa4ZZlRlFmHTfV+hkqwnbhrAIqxYWzv4V23VShtOnYIlBxz3jhO13cD+fcvMcVZwXB0b74/vZBnbOLI+r8rZtI/Zxl6vw6JJ6VTiHZCxDzrFafC8/Nkr1Nd+G5LbLCpdt2OxKBzLaiqp4XFJ0DRTklRzYdU5hHUjgSPRKbR+lI90iWLl9L1H74mvhr5VArEk9BtJXHgJ+PkJLpw5CiGCPQmNF9iD8gYE1AwvzS4qGK2ANc0hqAliqjm1ZiFMC1wI3e/zcCH4mINwoReSzwc1St7yWq+pND4ReUXJSdB+JNPh07tHTjFd6o1nezWTQK3ppbPM3r4nimmGPJ3XOVb1cpHOvCNejcpHeehTPj2EaQSUmhzXOFJdIgbyuTJeBjE36pufO3xH5c5u7Zhl8G3I7L3JG4dStUWBYmv9ISu6cUjFtZivdsyLwU1LhhfrUQ9+z6qYp7Ft/P1JGzdJQ4cnfhtL6axtWoeuHUc4tFq8moeG51tm4YbUUhUzdhc80kUxdOM+24kRulkCmSG2LPPaVg4uTGL89LpwwWxu0oL9xI+chcnXKUFe4aFXuafpGVnNHq+QKrCDJ/FlIPMXzir9Jb1s/e7xlD/LXbsqMUcik5MmnngZnCGboKpWkWCri52UhoDWK7i8Qh7mmYQ+2AMGDKbA42a16x/7d5JaQUgm6+5cD4ndc8aO3wEaMAPgr8lYj8AVU7vxR4q4j8RwBV/R8RaXRgPizzS8CjqZTM20TkD1X1fX1xMuCWYo9jjFdqCZ2Lqgpt+Vt3R979FVlS1qN951YXvB2559SjnmNTxIU3skdNscsSVwXWdk9JOXCI3CkR7Wr8c+XCjdiPbSMoF04BnLeKoMw4X1b5Hhe1crDPx5b0i5xl0VIARcZyacneEHshqAmnS6s1pfoDxPu1N3nUbtUfQLYMuBV03PzfzFSocyubz9WvNpUGplmYqcKoUrAkXw/x3bNP8GqG7FZXl7m4Z+e2yLxn47eA0oY7sn7qmkm5qARbHgHmWRfmRRaKHJmRs3FbLAoW5vnMYml+awVwgXG7IF9yNjd1aU7Ylyo1mRrFkmX16PzIVOBZWTpCP5tVhzOPpOCs2GeTrxTu2SqFMxTuuVYOZYfEj1A3G3H6kRq5SMOv8q/7bxYY7K2L0ptWFo2dbhZ2s4R2ZuYlNe/Y8IVnXnVreN5AtGi4NWcAx7poDDarNDLOm4ZVGktDpmXQ4uEjRgH8rfmz+APze5uIuEN4CHC1ql4DICIvp1IuvQpARDgS4Sjg1zfJK1r2AN/O68exu2YaWxUtsXsVZffB566i1HOr4p4XPFu9YTxdOPYpbc6auedigPRLzbo2dqQe5XumFp/4oRrFW+I/XyyMXz0rOLe04TI3sj+/tCP8muwL41YWQmncsL+FIIa8HYkfS03e5vB2dlwrAPd7XP1VRaXOLXdxTL0cq3t2v+dL8mMzSja/2XGJHBfGzQiwLJBjk2Fh3IrSPWtZ1n4d+5EHZ8/JEEOc5KbeFrl71qOF+9ULTLmdMW5nMoojU18XmLK9IKM4Y0xoF5j6PQOFebZ+xQWVe+Vv4p6B8owhnzOl+z0+U73b+aMq/6OjgguOqjJwCn6Rdc9RSMnCtKFFwObs2/Et8VvSP5sd10pBar+a+K0JqOTIKhRH+rh+XRO7OEK3ZJ4j9UDLI/isNdLtO6+SSzNcocPmIWeGFa/furjadRN1ysIqitLEt3Gsn13AN83ZzADq/m3djp35y5aBct7wha2PwpPDfz4/+HZxu4D+y1iYNXE3mqaka4EvawcSkcuAywAuvls/SWb0KwEfOTgNbZtCQxEE3FZFbiz5+0bfNkt/UbftZqEqbjHUDYjVs8t74erhL/WvLTB/pN0y40jpjcAbI3vTwD3lkFll4P36xA+QH5eIec7OV5HluKiJ/9i6LeF8RUi6tEphWSsApxQK1CqAMqAIzOhB8hw15COG7FkskMXCuB05Py3M89Io/2JRD/09qCEmq1cQcSYif2bhxhW5V36u7E04z9Smbi2l7Ky1hBB7bqbPrh4y2Vj4C9ht/yzwnCMN4q/8Mk9BZF6c5ju1ib73PQLhfKWQOUJWb23OrmmIqxBfEbiZh3tFHeSTWgJ1iiKkdH3zkS3/MZ4ZM3/FzAB2hVAr7K7pqb4I8wnKB9z/SI+HRmgGocXRhrlHB/xCWxw9bdy25fVO0VoLw+1nGz7mZHEmpWuImZeGawSmwS2ywtnd7YJfkYVKo7utLxMlzywJGZNDljmbc2FMPEWeUSxq0w9AuczQI0MuS8/Pkrw3O2ibcbJlPXuo3ZrmoCqNWhn45qGQm1UolhizpSKW0D0TUO3WMgX1wDWJTGpSdiYgqU02uSVscaad0py8Kxd4ph3zeySURy23hWcOOlLnZkf7emRNQQrGBGRNQYtFyeKoMM9m9L1YctbMAJwpKF9y1ky1zph2ckG2dLPg9s4faLZxO2u261iZlp1wZwN27SMpOKe16QfgiNre78/GzWt6o/0Cu3TSnh34bj6TtGcHMWiP7Kucmw0kZBZq8owlbs/04zXDNteU/gzA45mObZ98kHOOTQMrNHNWgj5MqQCuBe7u/X8xcN1QhBK4eYC8QxhbsKnSCC3MdN1CZF+Qdd00ayzYtN2crU6zxir+KjiSYnAK428TXLZ2f5wpC7fga81CyzxrLP6CMQG1F4HLzK0LWNmLInPrAYVVCmVWj0KN27Ks1wXwF4PbC76FdGz2lJ7y8GYe/uzCpdGy6Yu3CBx0G0N7DcB3y+r/a2VQmw3cGoBTHqC5NuJqXjrjtvUj03pB2P2W7tna5xeL0i0I+4vBC0fyZjE2L9xOHkv2Z/Jlg/irNAoWblePITWtFx2PjYbKM3X1VnqbGnzbP8DNcqax+AtwROH8611AZWddoFoYbs4ommcJmorKFK/nTy/GetsQr/StI0JrkTewR9/nl87W8xX45ZjWJhNvF5Al/YKsw3VtTKkA3gbcR0TuRbXQ/K3Ak4ciFAg3alPk8RO+AQUwsD20b0uoDRM6LzBE9v4WrvY5AX/fr195Q7MCt39Zs3pnhh3ll/4OCUMGZLUCsI0lqxWAHQH6u4DcdtDAziBVcW7+Nk+nDHy31i6ghjnKKgeVxnP165ku7Iv7Jipvxw+t+lffzSf4ENkPtZ2R7Z/1LiBDzuLFMX4i6j3XbraOxN8Gap7r33pYUm8D1caWUKhI0CkD5+btDPKUf9vtKCsaz5VfETwk5a89AdxcHlFklvjNOoMUZHLknsESu1UAddtsnw04kmVwW6fbJRfYLWThK4ChQ2TNOLHav8bQzsOhg1uhreV+uNDW0CEu8f19TrF90M0A6K7xtNGrAETkFxgYI6nq9w6mPAJVXYrIM4HXUY1/XqqqVw3FKcj45zL8SYIx22V4JhCwmQf294fCx1SKlbkK1yX70CndkMbvezfbKUrTAauTotbuYX7UIwgsmecctaefWX1raOhswNABr6V31sD/1YG4/qGv9mEvPy4BN9883w43djgsBJvccCiCxOifE6gPbnX9HekHDpGJ9+yHs4RtFX3wvIBoZ8TuHxjzD4TlA271jLF7IKsgq4e63gzTKm4b91jy4KGv9gGwSgF0R+/tcwChQ199J4bDV0kMKIAdHgQLHjSNOEBaubcXqcOHw9rnffrOL22yC+jtgzG3AFW9HLg8NnyhGR8vbr1SHjH3dowdGAud7G2EDSqD9uyhO7MIjQZ6T/S1dgG182vDdSx/kctTDlYe/1K49hUPCyk75OmfDq7lPw4qiva7+yeG/fRC4dpuPnySj703aNsIXvsQGC8NnQ5uk3k7jcwRdSCcR9zd6ybCxN5ONxelfRVEqVK3GX+HjD3EaGZmDRJ3s9I6P/90sEX7lHCVTv+Ivc8vtAAdcwdRCJsogjFuib3EsW+GUPuHBpb9HOHcvPMCfehVAKr6vwZjToClZvxzccte/1UvaVr12uVt3+kRih+6v6Pvvo9OhQdO/oXeI6d0Hbr+AErpjNLO9isavNCtVhp15wndIxR6x5BCaSP24rfxWxZXX/xbFUNmht77d0LKo03iA7to2vm2Dz81lciwScSVf6Os7A6TvE7PLbx6Gw/scsUI+XYukuu5Hyjkv2q4RpyJv4w21P7WuQAuFG7ogriYDSYxdwHdEfgh4AvwPgmpqo8aTX3LKMj5+HK1GcBwepsTyLYqcjBcYPbQ57/uLYGha6GhJpfCUxRlS1EAHWNhJurSGbpeup030LDF+DedhjB82VXcF5zGFMk2Lg4bu6VyKI/gPTeRMtUKXgPtubtrJxd1JkWXv7fbZWmVwoZXPw+/7/C7xYzsV0lvG4iddY6R8vCFjbGD0vDgK4SYReCXAb8NPAF4BvBvgH+MiLd1FCr8S3GLHecRN2rc5HrpsTRiKzpGllXuCF/1jv7BfEct6jXaZOLnOdp5t9C3p7pPPjbfWuF6hD24gB1Ko3Z0904FtgguI+UaG11v+0M0q+Q9F6z8RbyRPrPqR5y2ch008Fmq+qsi8n2q+ibgTSLypoh4W0ehGZ9aXrD9dDc4tLWKvTn6M2+R5BnzqblNP7Iyp89EboKpPxLjoybg1WQKETZ0CTt0hq1XlhCZRsXPT9zXvHaJVdvftj4Lu/YuIA/2K+wfE5EnUO3VvzhGuG2jRPhUsX0F0JvfFm3ImxDZLj4z18ljjXdd5522TcT7WPDdNraujALmt/i467fx00biq2Ib9bzaAHP1/GIUwE+IyIXAs4FfAG4LfP/KOW0BpQqfKUI3AW0j7d0vGHbz3DIZbolYtjbq3gM5H6ICOAmY+2cYN8G+2tQ2OGfTPh9zF9BrzOO/AJdslNuGKFW4eUcKYBNMTUL7zH8KRQnzMuEksJW1lykwdV/1sfUB4C5mACLyecAvA3dW1fuJyP2BJ6rqT6wu4mZQxF1dsA+cBNKZU4P3MVe59oHT/O6Hil19XjQGu+ShGDZ9MdVXwX4FQFWvFJHfAvavAFTc1caHhpPe6afsIHPDSRg4nCQcSt+bQs4YBXBLVX2rND+usNyRPINQ6ovKEqbDoXSoOSKV3X5w6AOSfbWTGAVwg4jcG2P1E5EnAR/bqVQ9UMRdSpYwfxx6JzxkJEWzPg697FZZnolRAP831X389xWRjwIfBJ6yjmCbQrW+iyQhYR0ceuc+aTjQteStYso2GbML6Brgq0XkVlR3AX4G+BbgwzuWLYjUgadB6qgnE6k/bReHNusdug76tlSj/7tRfQf4Deb/HwDeTXVFxF6h3k2UCQmnDYdGLicVJ0lpDs0AfhP4BPCXwNOAHwTOAF+nqu/avWhhpE6QcBpwkkjmNOHQ+GlIAXyuqn4RgIi8BLgBuIeq3rgXyULQ1DESEsZwaCSUEIdV7niKxZACsHcAoaqFiHxwUvJ3sqTGnZBwSNgFcSVsB0MK4AEi8knzLMAtzP8CqKredufSJRwEUgdPSKhxSIPUoS+Cze7ElZLIJiFhWzgkokrYDfZ3sc4ekRp2QkLCacM6vHdwCiCRe0JCQhuJF9bDwSmAhIQ5IBFOwknAgSkASR0vISEhYUuY5GIdEXmuiHxURN5l/h4/hRwJCQkJpxlTzgB+VlV/ZsL8ExISEk410tWaCQkJCacUorr/jfUi8lzgu4BPAm8Hnq2qn+gJexlwmfn3fsB79yDiNnAHqusz5o5DkRMOR9ZDkROSrLvAHOX8HFW9Y9txZwpARN4A3CXg9aPAFVQFpMDzgLuq6ndHpPl2VX3wVgXdEQ5F1kOREw5H1kORE5Ksu8ChyAk7XANQ1a+OCSciLwZesys5EhISEhLCmGoX0F29f7+ewzHrJCQkJJwYTLUL6KdF5IFUJqAPAU+PjPeiXQm0AxyKrIciJxyOrIciJyRZd4FDkXOaReCEhISEhOmRtoEmJCQknFIkBZCQkJBwSjFLBSAijxWRvxaRq0XkhwP+IiI/b/yvFJEvmamcTzHyXSkibxGRB0whp5FlUFYv3JeKSCEiT9qnfF7+o3KKyCPNFSJXicib9i2jJ8dY/V8oIq8WkXcbWZ86kZwvFZHrRSS42WIu/cnIMibrLPrUmJxeuEn70yhUdVZ/QA78LfC5VB+hfzfwBa0wjwf+mOrrZA8F/mqmcj4cuJ15ftwUcsbK6oV7I3A58KQ5yglcBLyP6vvUAHeaa5kCPwL8lHm+I/Bx4MwEsn4F8CXAe3v8J+9PK8g6lz41KKfXRibrTzF/ozMAEfmpGLd1ICIfEpH3mNHc243zQ4CrVfUaVT0PvBy4tBX1UuA3tMIVwEWtraX7wKicqvoWrU84XwFcvGcZLWLKFOBZwO8D1+9TOA8xcj4ZeIWq/h2Aqs5ZVgVuIyIC3JpKASz3Kyao6ptN3n2YQ38CxmWdS5+KKFOYvj+NIsYE9OiA2+O2KMMlqvpArU/O3Q34iOd/rXHzERNm11hVhu+hGmVNgVFZReRuVGcyXrhHudqIKdPPA24nIn8mIu8Qke/cm3RNxMj6i8DnA9cB7wG+T1XL/Yi3EubQn9bBlH1qEDPpT6PoPQcgIv8O+PfAvUXkSs/rNsD/u0OZQhf+t/eqxoTZNaJlEJFLqBrrI3YqUT9iZH0B8EOqWlQD1kkQI+cCeBDwVcAtgL8UkStU9W92LVwLMbI+BngX8Cjg3sDrReTPVfWTO5ZtVcyhP62EGfSpMbyA6fvTKHrPAYjIhcDtgP8G+AtcN6rq2NQnLnORDwKfoGpsv6KqLxKRhwHPVdXHmDCvoJpu/71ccOZBR3fp3GeUkJCQcLjYg344/6GP3qCBy+B6ZwCq+i8iciPwRar64R3J9eWqep2I3IlqdPQB4C3AfUTkXsBHqUZOj1HVqy6458V6lx971o5ESZgE8x0cJSTsDntu93/31B8OcvjgGoCxV75bRO6xC6FU9Trzez3wSuAhqroEngm8Dng/8DuqepWIPGMXMiRMBCGRf8LJg0T+zQQxdwHdFbhKRN4KfNo6quoTN8lYRG4FZKp6o3n+GuDHTdqXU22dclDVF15wz4t/eZM8EybGjBp+QsIgTklbjVEA/2VHed8ZeKVZIFkAv6Wqr91RXgn7xinpQAkzQ2p3K2FUAajqTk5aquo1wGQnYxO2gNTZEtZFajuzQMxBsIeKyNtE5FMict4ca57bNraEXeEA7JgJO0asXXuVv4RZIMYE9IvAtwK/CzwY+E7gPrsUKmGPSJ3x5CDVZYIPGT/KEfVBGFW9WkRyVS2AXxORt2wqW8KekEhhfkh1ctiIINZDQYwCuElEzgDvEpGfBj4G3Gq3YiWMIpHIfpDKeVqcILKdI2IUwHdQrRU8E/h+4O7AN+5SqFOPRDrrI5Xd9pDId2eYy+0QMbuAPmxmAPcEXgH8tbn9MGFdzKTyZ4NUHnFIhNyLuRDqoWFUAYjIE6hutPtbqq56LxF5uqrO8ha+WeC0NMbT8p4hnHIyToTLiWgDMSag51Nd2Xw1gIjcG/gjZnoN695wUjrASXkPHyegY47hRBPwSa6/mdVbjAK43pK/wTXM+AMHO8HMKq2DucsXwgnp5CeGiA+5Pg6wDmQm5R2jAK4SkcuB36G6tvmbgLeJyDcAqOordijf/jGHxjQHGdqYSYONwUGS8qGU78zLdi7E2sFMyy1GAZwF/gH4SvP/PwK3B76WSiEcrgLYZ6VM1QBm1CFmT8wzKqtezLAMZ0G6E5fLLNr2GvUQswvoqWsJM1fsoqJ2Xfl76mCzaMQwfyKeSzm1MDkR77lc9tJe91imUzSrqJPAB41tlOq2amZHjWmnHWFqUrGYKelaTE6+bUxQXofYDrct8tbbwY7b1clUAJvU6jpxN6ikrXaaXTSWyafWMyPWVTAjpXWIo+VNRd5K25m4b++6/cecA7B3AM0fsQW+asWsWAkrVfwmFbzFTr13op0ROe4bszG1wexIG9Zsi2u+xyp1sapc69Tzuv1w3XgxM4CrReT3gF9T1fetlcs+MFbYo/7jBThaobGVENkwNiLlLZHMrMgqYRh7UOIbTa5XlW+F8DHtNDb/2DY/lt64/3byWTWcjxgFcH+q66BfIiIZ8FLg5ao6j28CDBViyK+nkAYrY6hge+JFVcZIA1iZfDeZrq4d87Awe4OSxrLC7N/EIZqYtjAIG8qrL+5wnO2mNxZm6PWyNeUcQswuoBuBFwMvFpGvAP4P8LNmVvC81iGx/SFI7iG3ZsEEKy1UeIFwwUIOhgvk0ZdPfzLD+a6Y/mjUCTWAbpPLRgg05KuxpOvHWTmGH3kLhb2NNFaBaVux7+1LFyrfYJtuh1tpBtANO9Smw+Hj0oiO6z2HyHtVxZNFyjKkKNqIWgMAngA8lepCuOcDLwP+L6oPt39edG67hF84XgF0KtAvHM+vU5ANv1B+Ggo6kN76jatXhsi464bbFKsSa2z4IYXR5zcUZyjXdZRDHTkcd2WFN4vBfpucVwve9NOwUmg7qHTaqio9iqEZW0Q75Rw7Yh/rlzLQ97O+cBGE7ofpSyfkH5IjFCaEGBPQ/wf8KfDfVdX/EMzvmRnBtLBvHiL9ANlL0M1Pr1u5zTj9jcWGC7k14/QrqHVGFz5W0f5jefhYhwjLgThDuYXyCruF/MPk0iYDP4z60miISJpujaTG8goNdEPlEiqQofLbhlJQ2dms0aUanDWL91j713H8oIFZgSkXabv7+QcGcKp+foNqP2pwFhrhh0g/C/BBH8G7dAJph/4PxfX9s5HRw6ACMKP/X1fVHw/5q+r3Dqa+KwhxxC9eAflkHyjkOlxN4kOEPuYXMzJoVl5AloG4ff5t7MNQ0Jd7SAG0yc8P4/uVAcIOxW3nHQqnWqdtf5sjRHHhXNvxwnfeT6WOG1JC2gzr0m65NcN18+j163Xbr1lomOQDEQIDrRDprzeok5ZbHXa8n7fEFAn00a7SUnpes0X8QYIX7ZB9NhLO97PEPqoARhT8oAJQ1UJELgGCCmBOkEy7o3zpNpxMtFH5Nry04obcMr/SGhXUypfuiCDzlEKookJ5uLQIh2vn1YgTctuSLaEMNPsQ2Vu3PpK3aQ2F8/2s9KVHzmWL2ItSOnGVrgIoS38UUXdyLZvhFGoSLz3ZnAKwfr6b5xdyo+0G0lYKSpfkVbpc26NEogb2fWEi9EhQ1wRIvBEuQN7+r3bMpl46fhqtuBpww+u/gwPCrO6XvnKwMw9/oJdlpXn2/Iy/bWuZJ35ppM4lng/a/TtE9j5HuF+UTMpu3E0UgMFbROQXgd8GPm0dVfWdEXF3h/YI3B+xZ3WYzKtUgCyrG0ZdoV1izUSdv0/wedYsZAEWLp1uwedepaxaabbiK7fSpOcRP22ZS+cXCufDD7suSq2be1spFI3RfGZ+xYUrHElnNclTE/eyzJpx8d1sGplzs4S9LDOXdmH8ilIo7bOTK6P0CZ3wTEFLaqVgfiml+qsEM4mIexbfz844nB8urqui0uNOG85XAH4T1lbcgKKQgJufzijGyJ02sde/HcUgWru5EZUXLvNE6xA7kGkzPwH1+rcLb9Oxbpn3upl6cbsKwHGJ52efXbZZSVnmVXLWT7Tx7LJrFYFq3Uz8funCB0b7Cym7bpZ7enhj0XErg/n5iFEADze//ixAgUdFxB2EiDwW+DkgB16iqj8ZF9F/9kbsGe65dquJv/otO5WWe2SfZ3WBWrK3hZhnZbCCbMXUFVDWleWFt5V25FXUUVadsbMkvcgKcrqVeySFkaVuBP6z9XOKgroRWITduopg1XUEf+ReeM3fknfhiD1rPAMca07hPdv07LP9XZY5xybc+bJqtsdlztK6Fbnzs0rh2HTY4yLn2Lgtbbgl2K7aWA+wyqgwpF8ILM07GTcpBMzRSHFuIMsW2ReQef5Vwi1lYOO23bxwtVJQxNrGGkqBrlsbK4z2VQJ+7RG9R+LqRkheOEfI4pGz5xdya8fNcEqjdtM6rslXfQVgFUbWdSNrPldx1cu39nPKwHJFJs5Ny1phlGWTXzRT1PKGKaqihNzmZ1+4RxGEZgo+8UPFEW0eCnGJ79aHmG2gl4yFWQdmfeGXgEcD11JdMf2HsYfNpN1IhYYygKqCfOIHyPNaK/qk3yb7RVY6N1vIeVZ6ZF+4cNZ/kRUuDetmK2AhhSNxm8aRFI7Erd+RFy7z/JwCwLotnaxHsjR+6qVXu2XYd6uVQ60oum4WuccaeeTwsfBYo032BeIUhHU7rzkl9TPAsS4c8d+sR8ZvwbnyyLhVzfam4gLOGWVwU3mmClcu+ExRhfv0snK7uTjiM8vK7dyyCi+inDuunks3ws9r4j82PfY4Q47N9P68IfNjIauK17llx9Rux3hu2nRbem4mvBTquZlyLiEr1PkDSOkpADfzUG824NVRu7oaZO4pvNboXHPxCLEmWM2l4VaRqY1Th9cWOZe5R955Tea1G86tzNvpNZ9tGg0FYWXphKMjS5VeSwHkGlAU6r2nccvVUy6mPtQ3sfh9xyToKYJSbRxxv2Ojc2iagCy/+KP9EOf4A82xmX7UXUDms5BfSHU1NAB9C8Mr4CHA1ap6jcnj5cClwFqnjUW8Nu5Pqaz29KZ0buSf1SN7n/itmyNvTxEsWm5nsqUb0Vu/Iy9cg9iN2wWGnI+yZcPf/p6x/s5t6QjYEvsZKdyzrxzs8xlvFnFEU7kdoeR1PzbpiXvOTUlmIu7ZIvNG+PnI3tTCEVIlU0lJYd7DdogC5dg8G47kWOFmpwyq/G7WBZ9WQ+hGKXy6vICbygsAuLG4hfk9y6eKyu1TC+O3PMuNmXGTC0z+wrIwMwUzVlPMiB/AKIDs5ozsXOWWnzPvfbOQma9iL242bue09j9n2tK5ksw85+dMHZwryM5X5SHVNKT6PTbPSzNiO16iS6MhCuNWlO5ZbdmWwx2crHoPEaltmblhyyx3bmLCsVjUw1XjposcFqaMFrWffXa/eYYuLPEbRX8kNfGb33IhjrBL54bnhglPJ5zv5v+WnhLqhtParayVQZVwHU58pdSaVan37KN06XT7gl0rKKmJ31abqjhzVekphZjZd986YdfCUPNBH2LOAbwQuCVwCfAS4EnAW0elHMfdgI94/18LfFkg/8uAywDyz7ooIF9/gW1y8q+vItr2+XCYMmjGiV2EzQJmGTd6d2af0imFelZQeGYeS/Zlg/irNODIpWvInpr4j0xPCJG975aNrBYunCncvndOaWQ9xp+uNqfMVU9rNtyCgrNGRViT0RkpOPaUJBilamYIWVnXVWf9xVtwayCwGOtMzWXtZ2fWbnRe0BnFZ8dKft7UlyH97HxRE//N1fvIca0AOK7cdLkEowD0uFYENfFbJimdScJCGlvKDDlngljiN2Ulee6UgRrCFtWaTQ3pC7XCEbXh/QGXnSlobRqz5byE0tn+7UymjuPMtqXn75t9nFnLvm9gB5ZXR6FFdA24OdH9cAxAe56HogyYe7aFmFnEEKLWAFT1/iJypar+FxF5Ptv5CEyIPTpvo6ovAl4EcME9Lw74205Z2/KE2knbmlfULQjW07bMTdeysiY4AsEojNkgs6NaoTSkXHomj6WdKZh79JaScZxVnceN7DUPmoDOyZF5rjpqJtqZFfijfd/s054p+KadM5Z0xVceNSHmNHtRxjLYwIbMQQ0TUGDHj/U/9sxCx57pBypTkB2V32xMO8eaN0b+AOfKI2f6uakwv+UZzwRUhbtpecRN1hxkTEHni9ytB1gIIHlVVm4ke0HpkWgVLjuC8oxRiGc9E5AxFeXnrQLInOknP66VQm0qqhWGLEv3DCDL0pl73G+hbsRvFU9lAgoNTc2vvyLZGhGp1KNzb8tbx9yjmUDAzY3ovelkuWiFy5vmoMrPmw2EzD3uVzqmojKv38mPSyecdtzUWwNomIzsDMFfH2itAZDX64n+b22zr4uxvYXUH2yEtoqHdgUOIbSZIhOtN0zYdDSj2MIi8GfM700i8tnAPwH3ipJ0GNcCd/f+vxi4buVUGtrd2nLtzpHQ9DhzU6/MKYeys/sjV+HYEIRvFqoXfxfOr22PW8dGV9v9tfFcxa3JPrTgWyuFrm3fb1TO3i9lZ5YRJvq4nUL+wm9j9w92503mwrkdPNQLv/6CsP21cayN/1hz5+8vAp8rm27ni7z2dwvDuVsktm7Hy7yzhVREyRZmWm6UpWZZRSbUpC9LqRd/zcJvVtQLvW4xuPDs/GXt5sLZ39J/rkfJoYVhNyL2FoaDi7+rDgztOEpCbtLxbyzaOiZr+WPIv+OGZ3evfzu2fQm4+SQeSs/PKxQuRPbezqHqV73F3/q3s8nE21AS2hnU3N1Dw20V+LvjnJvROEurBTNvYFbavEq3vtWHGAXwGhG5CPjvwDupmtZLYoUfwNuA+4jIvYCPUl049+SomFqP6OtTfuLN9aztIXO2ZjuFLUttrNhDVbCFNRdY80dWukpbets8a1L2Rs6tHT8ZAbeexZyYLZ9NE0bt1zYpZVJ2GphP4qscEInB0N7/5m6gemYUsw30uKzJ2e7yqbaG5h03u9XTuhVlvTXU3w5amI5gR/2lSrdzeCO1zPSMUsp69Gvtx4W4MwGWccpC6lG3JfvGTh7PzRa9I3Hv2d8i2iJ734RRmzwCtueAqaMRZwBj+/uH/FW8sKNuTXJubA0NxPW3fvpxnF+bxMUjeanj4qXjwnubRtyvHdH7ZO42j9Sjfn8reeVXb4jOvU0moVlBG6W3BuBva7YD0KWZFWfa7cdlKY4PCjfozDY7CVwJos8zj78vIq8Bzqrqv4zFi0h3KSLPBF5HZf59qapeFZ9Ay1aL1vu0XeXWh4fciE20oQwqt8C0LctqN1vwgalcJt1K7TvlF9riFVpTCJ1JwEvHPQcUQCjckNsmGFIA1bOdiUnHvxmuORIPHQ7zD335B7yK1tmAKpxJ1/NznOvODTSfgcbpUYss94jEtq8cb7YpdRouE9+t1U6V4UNf/j5/39+gY+tux2ljrMrb/sMDxqB/8LBXSGkEZxdevHY4od7XH0zD85Nxt7HDoo2+3TLtZAES7zs/FD7p318R/iEy1y9svoob0Ph52MFLPXCsB392MBSz7hi7C+jhVBfBLcz/qOpvxMQdgqpeTnWh3AZpVL9SLVVVbm4GoK4GnSKQeqbg2kegwvtOB+Pi+OGa6QEdsg9d+xB7+VPDr8d9yK2NsXt/htB3J9CYMnDxB8K0T/A23boyaEiJ9Pj7uy/sb9eNuh34Qrr68tKybcjaklU6JN53/UPHZO+Xk5Opxz+EdYh/U4zMFMJxfE0WiNtulxIKF7inx1MewVsA/Li0w4X6uRfF67/B2wJMuKwRrh4cNvLCH+R4o3e/fu2gxVcsbjxbDxa7M4m8u9EhohHE7AL6TeDewLuot2YosLECWBuN97IdX+uRUGj66YL3KIVAx+8QZaDxSY9b+DkyHF2EiD2GyDch+3UwdGlcrPII8llAKTT9+/MJKYIxWX3YUOoTRDtu4AY46csr+ILD7zQUrhtmPMjaWIXkrVNQWQTI3oXvVxRjFzc20mjlEXtxY99ArxOOYdhBS4Z0VtMy6Lgp3b4Qmv33DRyHrAl9iJkBPBj4AtVgc5weTqqaxBuzAsfy1k/CSsEqkkaDaLe+cIMLFvJA4xuKt9oV0XFVMvVXvWJbzhghDyuXzcOv8i2B9r01QaXQiDBSCAH/kDQbXU+9BUQPKAbEjP0mRzBYQzkM95+2rH1+Q+GGULOG5xaaCXuDUzsrKLR71XVBd6DXayVo5bGuNSBGAbwXuAvwsYiw00GhXSyha36lthQ1lYJFYKBWV4I3khxSFD4iG3adVk+l7YHsdz1b2BZ5rX6ffly+G8s3UH4yNl6MzLp9N+k6w7KtDQhi22RMmBXa/SqDpL7wQ1DtxlGPsJt8UbnZDQcCnfbmE7Ed9Yt01YeIum3S2/j4S0x/jlEAdwDeJyJvBc5ZR1V9YkTc3aHLyV23oA26Rym004BhBeEnY72CBd6NGwxl/cZIaEWS3vZ4cV8mpSlGujt9s6Fy2+Bd1/n63DawisTxs4YYk+YKGQ+geS14cK41EKfrVssXiBcUYHhIsA6xr9M3YxTAc1dOdZ8YUgRtdxhXCoGwvQoCGo12kLT8mUcoKxtstBLDs5WxdDuprNmRdjlK3gsmNqFsHXsqz72UmtsB0/9O/oxniIhbsQbSk2BfGCL20KwrODuofaP7m01n6GNK0abfiDAx20DfFJXb1AgZ5Kx7GxFKoQrXXFMI5xtZuYGZRTBYrywjaUZgrBGvjVXY4YQR8NTrK2uV5xpKIyZGnyQxg4YGqcWa7CI7gJ9/mOxH5GH8HQbNwMHwvUOziLBxecWI1KsAROQvVPURInIj3TG1qupto6TYJ0KzgaFwPoIEG1GC0v326GiUdcZTAzOP6CTGZijrYi079Dz3FAwikjymxmjH34YSjjV1xGIdmRrEHqsMAo6hd4ndIeZMt1HB/Zgr9IGh2cCq+TbRqwBU9RHm9zabZTEBhkxAMXFWiRvbcLc0+nYVvu2RtLJX80zfwvsc4TrqDMl+30ppV+2v2sOxvuBDhD0qaTTZhxRFVNRNso3GOsog5hzA7QPON6rqccB9foidFQzFHUK0gllBgEgb6CYINpZ9mmcaU/79ZTuKCHvwVAjfXrpfGXaaXcfuvqVkhzxXUDrrtoPeWFtsV7KGJQLiFoHfSXVp2yeo3uUi4GMicj3wNFV9x+rZToB1ZgWrphvCWopnTQFXaszrZTEqwi4U4jYxVkZTK6OBYplSEU0xC9Jd9dlGJtuzp/dmsX7U6D69btuIUQCvBV6pqq8DEJGvAR4L/A7wPwnc4T97rGvu2VZeIWwj/20RxEZT8e2IMISNRoZzVTwWUysgmOcsyJTLVGtHve1618Wyiy3CHqJOAqvqM5w8qn8iIv9VVf+jiPm80knAPpVCbP5D2KVsWzdObrfT7nvhdSumiKkJdJU6mIMSsmgV2+SKqI1WWc1pc0NMP4lRAB8XkR8CXm7+/xbgE+abvnEXxh8qplYKQ1innU0l+6467d4Ope0lmw62usV0TsR5qMoohLkrqBHEKIAnA88BXkX1un9h3HLgm3cm2Vwx1CDnXvfb6Exzesd9dLYJR3RTbzHd2RmHuZLkOnU9dwU1gpiDYDcAz+rxvnq74hw4drEgPDds2uAPrQz2RVYzMh1YTK2AfOzlwN1cFVMsdnEVhIjcEfhB4AuBs9ZdVR+1cm6nHdvcVnqo2DapnJTymvvi9MSYkzIKYfIT4bBWG4oxAb0M+G3gXwPPAP4N8I8r55QQh9Mwi9gmprrS4qRg6lHvgSigMcxdQfUhRgF8lqr+qoh8n7kX6E0ichj3A51EJAWxO+yiE6f6GMbUCghOjBJaBzEKwJ74/ZiIPAG4Drh4dyIlbIR9njtIGEcyec0fc1BCfdixcopRAD8hIhcCzwZ+Abgt8P07lSph90jrEYeJXfFBqut5YsfKKWYX0GvM478Al+xUmoR5YarTjwn7R1IspxIxu4DuRbUN9J5++Mm/CJYwHQ75LETCfpFMYLNGjAnoVcCvAq/mpJ/8Tdgccz49nXD42JZCSW0SiFMAN6vqz28zUxF5LvA06u2kP6Kql28zj4QZYR+3OiYkrIKTdip+TcQogJ8TkecAf0Lzo/Dv3DDvn1XVn9kwjYRDw3pft0tImB8O6T6uHsQogC8CvgN4FLUJSM3/CQkJCQmxmNk2bdGRI2wi8gHg/qp6fmuZViag7wI+CbwdeLaqfqIn7GXAZebf+wHv3ZYcO8YdgBumFiIChyInHI6shyInJFl3gTnK+Tmqese2Y4wC+G3gWap6/Sq5icgbgLsEvH4UuIKqgBR4HnBXVf3uiDTfrqoPXkWOqXAosh6KnHA4sh6KnJBk3QUORU6IMwHdGfiAiLyN5hrA4DZQVf3qGAFE5MXAa0YDJiQkJCRsFTEK4DnbzlRE7qqqHzP/fj2HY9ZJSEhIODGIOQm8i4vfflpEHkhlAvoQ8PTIeC/agSy7wqHIeihywuHIeihyQpJ1FzgUOfvXAETkRvqP9aiq3naXgiUkJCQk7Baji8AJCQkJCScT2dQCJCQkJCRMg1kqABF5rIj8tYhcLSI/HPAXEfl543+liHzJTOV8ipHvShF5i4g8YAo5jSyDsnrhvlREChF50j7l8/IflVNEHiki7xKRq6b8OFFE/V8oIq8WkXcbWZ86kZwvFZHrRSS42WIu/cnIMibrLPrUmJxeuEn70yhUdVZ/QA78LfC5wBng3cAXtMI8HvhjqvWIhwJ/NVM5Hw7czjw/bgo5Y2X1wr0RuBx40hzlBC4C3gfcw/x/p7mWKfAjwE+Z5zsCHwfOTCDrVwBfAry3x3/y/rSCrHPpU4Nyem1ksv4U8zfpDEBEPiQi7zGjubcb54cAV6vqNVqdPn45cGkr6qXAb2iFK4CLROSuexQ9Sk5VfYvWJ5yvYLovqcWUKVTXfv8+sNKhvy0iRs4nA69Q1b8D0BUPKG4RMbIqcBsREeDWVApguV8xQVXfbPLuwxz6EzAu61z6VESZwvT9aRRzMAFdoqoP1Prk3N2Aj3j+1xo3HzFhdo1VZfgeqlHWFBiVVUTuRnUm44V7lKuNmDL9POB2IvJnIvIOEfnOvUnXRIysvwh8PtVnVN8DfJ+qzvFK9Tn0p3UwZZ8axEz60yhiDoLtG6FrkGLukNz3dqZoGUTkEqrG+oidStSPGFlfAPyQqhbVgHUSxMi5AB4EfBVwC+AvReQKVf2bXQvXQoysjwHeRXVx4r2B14vIn6vqJ3cs26qYQ39aCTPoU2N4AdP3p1FMug1URD4IfIKqsf2Kqr5IRB4GPFdVH2PCvIJquv33OfmDbkk6fpCQkJCwCm7kEzdo4DK4qWcAX66q14nInahGRx8A3gLcx3yK8qNUI6fHqOpVt5Xb65fJV00pb0JCwtww4xH2XPCG8nc/HHKfdA1AVa8zv9cDrwQeoqpL4JnA64D3A7+jqleJyDOmkzQhIaEBkfn8JayNyWYAInIrIFPVG83z1wA/DqDV5yEbn4hU1RfeVm7/y/uXNOFUIxFMwgnGlCagOwOvNAskC+C3VPW1o7FSh0xISEjYCiZTAKp6DTDZydiEhIQ9Q+aw6zzBx9SLwKsjNaKEhISEreCwFICAZMkElJCQkLASirDzYSmASgNMLURCQkLCicCoAhCRWwLPprp862kich/gX6nqNN/xTTOAhISEhK0gZgbwa8A7gIeZ/68FfpcJPuQuwJyPVSckJCQcEmIUwL1V9VtE5NsAVPUzMhULi0CeT5L1zjDhVRwJCQmnGzEK4LyI3AJzOZSI3Bs4t1OphpClNYCEhISEbSBGATwHeC1wdxF5GfDlwHftUqheiCB5UgAJCQkJ28CoAlDV14vIO6m+FCRUd5rfsHPJQhBgcWAbl/aNZFJKSEiIRC+bBr4L+jHzew8RuYeqvnN3YvXhBK4BJCQkJEyEoeH0883vWeDBVN88FeD+wF8xxYcYRJCjo71nm0bVCQkJJxG9CkBVLwEQkZcDl6nqe8z/9wN+YD/itSACizQD2AqSUks4VKSt4FtDjEH9vpb8AVT1vSLywN2JNAABPUprACvhwIleysOWf47QQz9MmRTA1hDDpu8XkZcA/5tqK+i3U32oZf8QgTMTmID2iSkJe6q8B0he04RvP5hKKUxN5lPnPzFiFMBTgX8HfJ/5/83AJB9m0Uwoz85kBlDuLyvZNjFvK71ty7XHMu3D1st6ADoH8tn2ruptvtOWy2cv5X1gu9RjtoHeDPys+ZsWmVCeHZgBTG3u2EL2IQIaTXYNM4msGmWdsl23PjZUBJuQ+KEZnDYmtXUJa5V8V5RR13mlNWcwG5XfvvX3DhRYzGVwHyTQL1T1c7cuzQhUhOLsFm0CW+rtKxPOisFH0x8jzAj5RhXCUBojccflb/n3VHGs0nLBdjUgWCXdXY06t02qIwQaRZRjQUbSGJRxLP9IRRZN+CtW21ZmFxNMCGPsKQ/2ns8C3wTcfjfijCCD5S22bxRea8QYG2WUHIfi9nsG4/XMBIJhA27BcogN52QIJBEI35ApdwH709VwUa4qSwgrz4Z2QOrRo96IkW4vGYWcB94lKFMofIB8gzIEnAaJMxg+EK6nTNZVKKN1sal/b7771wAxJqB/ajm9QET+Avix3Yg0IEsmLG+5RSObI4+4gh8fJa8Yd5RcpT9ciNy8cI10IvJuyOfF7chdqpOrGcf8evpZ2uXbk0dN1NLxEz+JoHmsW3cuX58YAsoxivh3MYsIdPS2LLFEN0q0AX8d8s/8cJ5fK1ifX4zcffl34krPc4R8ABqiilglFKt8BsJHxw2mtx9lEGMC8k8EZ1QzgtvsTKIBaAbHt1x1bjbkOd4RV0kvfgQ9lK50/RvE7rs3Azb9JOxu0uvI6udRdgnbhVcJK4+ALFqG5OrmZ6vB5VF6pGAC9pVZ5z08ondxGsowpHh8uUN1uEUlECTkQDv0ycv6l+qebfUK2iVYZJjA/PyylgwNMg+52XgSIGzPLUDeIcJukHRIyXTcwuFDsozFGfMLptsbbrWZzGgeQ9iibogxAT3fe14CHwS+eXsixGMtBWAQbeZZx2TjkXdcuP68gmTpuzfc2q05TOyDcTWUr0e62g1fKwPjV+JasfVT9d3qdO1zQzlYf+uW6bD5xuZblF2S1zqu+G5OuXh+vj80Zwltv/bzqmiQnzR+BWoS9/wcqWRevFZcMmqF6fwUdZcmatPPR1YTlzsbIE2SB/N/S/FUH+cwOXiKwhG67+eea7+QMmorjxA5q5dvHW8NRTE2axmajawSry/cSJwQdmEiilEA36Oq1/gOInKvrUsSAc1gect1Y7fIKCrDgdSG0ok294T9h2YA4VG57xYavbfi9ioF6+aVlVMARimULf9Weu7gVtBNnOmqEa7tVmrt5pNz0QyneQZF2ZSlBClN5MLKXLpwNdmXiHVz4UsoqkzUuRV1HONHqfW6RhnQVObKchGpid3eYZXlzk3s1eZ5DpawjZvmWe1vO36eoTbt3BAymTO7uarPM2+EbeNKh+w1oyZ5361F7JrTSc8n4qBbkOzr8INEHCL7sXC03bqzjc5zT9yOe0/cXrehNCLirZTOhohRAL8HtC+G+z3gQdsXZxiaQXF237n6AsQHjbMvb5CWQlCpBc09cW6d2UBDUQzMAPy4pXTiuvUKn+wbCqUZTnwF4Lu1FICU3rMje63dlp5CsYrC/aojeRuXZeEIvaEc2sqjKIZnlG50nnVH7HnWIHkXzlxxonkdVxcmnCXYPKsJe1ETd+3vu9GK2yJ589t280fxtcLAI2/r13ULkr0/Og8Qe2hGYREe7XfdViHubZpbVkprhyS+CYZuA70v8IXAhSLyDZ7Xbal2A20MEXks8HNUY5iXqOpPDkbIoLjFFu2xW8AutfO2sfL6xuisZUABjSqZrqIIK5RAep75yPevwmvArZuelIF8Sx2WpfGeEe2wx7wRGsH6JpjazTwHyNkn307czMsjkrBDtvhocg6OzuuC3PloemZYeVfZhBiaAfwr4F8DFwFf67nfCDxt04xFJAd+CXg01XeG3yYif6iq7+uLowLlhjdBRDekUC32xN2KDdHPbyhuo7PpoH9XppHwMel5ckqns4fzl1H/sLjbQpCrezqpDjWQTTr2QLKyQlubuqyC5RMYFGjALRw+EGbMfzD/6kfG0gv61XG2soY3EK6d3xB2qVCGbgP9A+APRORhqvqXO8j7IcDVdn3B3Dp6KdCrAAA07y+NeHKPcwsTeyv/EFmGSLrHv5NuY2TV71atBzYJW0S95zp85sLVcbNM6zhU/4s03XJRsqwaMjtTclY2/AEWWUnuwlm/koVxW5hh9yIrOm5HWdF4rvwKjoy9x4Y/8tzcb7Z0z7npZUey9MItATgjBZnJ44zxyygbz1Zmm44Nn1O7WWTBPbhhlK29iAVCYYbl1q9AKMzzsTHoF5pRmnDnzYVIx7qgMI3yWBfmN3fPBVK7lbV/Ha56XpaZc1ta/9L4adZ4tuGXATerDKxbUWYubVtiRZlRlPZ9TZmUGaWJWzqFIZRl/VyFk1qReOG0rWTUU0xe+HZcfwY66OY5hXatBZWLbxXtKIpuh5dQwB7lEcP/6yqJIRPQD6rqTwNPth+Eb8il+r3rZelwN+Aj3v/XAl8WkOMy4DKA/Ha3W38qODgy7jr1js7bcTYZnfcph3Zcz60x6vaJH6qtg/aZ2k9aZF8pgLKRXpaVjtAdwWfqiN0nfUvKuUfm7tm5dcl+kRUcec8AF3gkbsOfzY67ZC8FF2THDbczsuSsHLfCLR2xWwVwJAVnaKdXcuSI3ygPFGeCN+WXA7kpJJ/K8wijbuF1XasyClW7lu25wbFJzxE8GcdaEzXAeXKP0Kuue95TANbvZj3ifMvtXHnUUAYAN3tuvlI4Z5SHJf3jLKwArPLITdxCMvKsVgZgmq4pX0vShUBhyD7T+r1FauIHEPGVQqvQbOJgTAMmj9oRMc/q6sGvsyE3agIITnhMP1EZ5pWBdFUCtCJ02b5HvE4Oa/LikAnI3vj59vWSHkVI5M4rquqLgBcBXHD3u6sUK75prGqU/n/HD9uE0msnGDKDeA2ob9Zg/VoDHD9NDYTrVS7Gr/NKAYWChM08vnJppydtRdVw64Zr+HviZC1l1H6uwnTzCCGLbQMj6ewKg2anFsqBsH467XNvvp+2Rt/gDWQD4VT9Z9+/G45WuIYYjXCt/Ly4jfQ7I/tueo1RfNuv4dYNNzSa783Pxe1mMWhuWiXuELa4IDJkAnq1+f1fW8utiWuBu3v/XwxcNxRBFMwgsCdAtONqWCOJbl2OJLKOmNtrB02MNURt/ISjBjrC6GJxYPtpaDF2cIFW6S4CB9z6dim1t70Gt8m2368PDaUecvNGmZbLYnfZhMJJIFy0m3bSG5O/szBMF6OHbKdqw1uIs5bo2x5fbJjekAno1UPJq+oTN8uatwH3MWcKPgp8K/DkwRgl5Oe20GJ2PMiLHkRGEu1o2rEjjRCZDRF1X/go0tX6OXDoK7g11G7f1NA2UBrbP2349tZQAltDs6LeGkpou6jdBlqWjWfnFzowFjoo5gpL6t/QAS9vqydQ/e8/Q2N7J55baZ89N6cMBrd8SmdXkTYOglnZCe8+CiqjppmkeYagDtdxoydcn1/Lf9X99Wvtx4+kmZ3vTtpx+kMmoJ/ZZcaquhSRZwKvozK1vlRVrxqKIwqLz+xKoM3DxR4OW3XquNKBsU5cj4jXIfahUbcj7G7cautli7DVI2obrvBG2z6xO/8uYa9F4v5hLxtuWXTc3EGw5bL2s272IFhRoO5sQGl+FHFkbw9rZe4AmNiDYHmOLBa1v3HDHQozimDRdSPLNlIenfMCUisK/0oI9xEeR+wSUADiTJWD5wACh77GzgHEKorQVRUrK4+BGU0w/CrhBsJvFG6LGDIBvck+i8gZ4L5U3favVfX8NjJX1cuBy2PDSwmLm7aRc1uQ9f3Xuv9nVf+RkXu8W02+MXEqYteOXzdcd7RfhWsReyi90KlfP26h3bj+iWEnSzOvKmA9IhbLcFIzmVpSDlwZUSuquoCaF+xFjBhC1z/4cvmzg84VD103zYWOTaXEFaa41eVaGam3oOpIyikA8QjTmwn4J3tN/iFiD5uoWrMCb21pKL2mv/d/28zUkNlvvK00Wmn3pTe+AaTHfTD+AJOPkfwWlUDM7CTmMrgnAC8E/pZKvHuJyNNV9Y83FXBVVAogdqgegcikYkf2g+GD4XTEvz9+yK/ppt2wKygI5xYy7bRl0sA7e/f+4M8UOjJrN+2RXZaug+XSIX4R8dYhahJXe37ES7t7b5J6jwOVvsoHa4Yur+1TELRIJHQxXOiQmSPukZ5vywr17hkyP0p9T5NLRgN5aEBmvClg4F2cwggRtwTcCLj5cQNl4IUfHO2H4oVkXiW+CzdOLKPkvKfZQOxlcJeo6tUAInJv4I+AvSsASlh8RmPKd3VsRRn0e8YrhdhwYeUxdk30UNxhJRNQKKGbN0MyDSi6hp8v+5Bt3UvEbfuzppHQddr+lpi8lreTQ/SFgSs0wOgPkHTDBUki5kroBqmFWEtrP/8CPqgKJUCcruSd8qBT6Q15vVmEU8mDhN1VClW4gXcz6fbdKLptEh828/R7bnQ+aQWssx4RowCut+RvcA1w/epZbQ4plaObdvvh2BXO+FQY4YLBO2P6vAaihBVE1zEYLvZO/DGFEwgX9A8qo1gZViu3VT/qE+osEtsDd3Acd5PFRHcDqzNvNVL2wrUjes92ScMnbBdOvCjGL/DdgMYrBAl0A0IOKr5AuEYe2h9uIP9Bd/pmWGEFFoPgdwv2hBgFcJWIXA78DtVbfhPVtQ3fAKCqr9ihfA1ICYvPbFEBrPEt3TZWmo1EzzIiAvYEiYrbV4SrzmAi426kJCPraK1Z4S4+9rImuuTcwyLFgMxjhdAh0aFZVk846110nfpkHjhTNa5MB8hx1NQVQcTb/kTkxruC1vy28bqIUQBngX8AvtL8/49Un4T8WqqWsUcFoOQ3LfeVXViGbXDGOsSzYpSVP3O5C6JdVYZNCHm3E8MGYst2s/vbVy8LGRtJFgG3loyjEq/4TtGHZVcgvuCMYwjrVMMWZnk73yK6BcR8EvKp+xAkCiXk50KteAVMOerbdt5bIL21vofcxrbea591sweFsVH/X8cssE7X2KY5K5DWKqnHKszBUNs2p+z69r2J84/ZBXQv4FnAPf3wWzgItjJEFTk3dBR4y9jjqHIMWyHqWOxbSc7IFDMbbDjOWRn7JroNlcUqmOJj66OY0O7vI8YE9CrgV4FXMzUlqiLndmwCOolkdGjvdGjyHirmSIxjWEPmA3zLGlPPAICbVfXndypFLEpFzu9xBnBSkAg14STjEBXZTBCjAH5ORJ4D/Alwzjqq6jt3JlUfVOF42kXghImQlNj+kAj11CBGAXwR8B3Ao6hNQGr+3zMUlkkBJFQYPK2bsHNIUhQHjxgF8PXA527r/p+NoFpf0pWQkDApkvo9fMQogHdTfRd4ktO/DSj1zY0JCZtAUzs61Rg9NHE6EKMA7gx8QETeRr0GoKp66e7E6oO6q3kTEkJIZqGEOKQBAMQpgOd4zwI8Auh8I3gvUOo72RMSEhISNkLMSeA3icgDqb7W9c3AB6muh94/NM0AEhISEraFoU9Cfh7VZxq/Dfgn4LcBUdVL9iRbB0qa4ickJCRsC0MzgA8Afw58rfctgO/fi1RD2MINngkJCQkJwwrgG6lmAH8qIq8FXs7kp6o17d5ISEhI2BKGvgn8SuCVInIr4OuA7wfuLCK/DLxSVf9kPyL6QoGmGUBCQkLCVhCzCPxp4GXAy0Tk9lQfhPlhqqsh9o80A0g4KUh70RMmRsw2UAdV/TjwK+YvISFhE6TBTMLEWEkBzAJpF1BCuoMmIWErmEQBiMhzgadRfV4S4EdU9fIpZEk4QKRBwHpIijOhhSlnAD+rqj8zYf4JCacLh6Y4k8LaOQ7PBJSQkHA6sEuFlZQLUJ3s3X+mlQnou4BPAm8Hnq2qn+gJexlwmfn3fsB79yDiNnAH4IaphYjAocgJhyProcgJSdZdYI5yfo6q3rHtuDMFICJvAO4S8PpR4AqqAlLgecBdVfW7I9J8u6o+eKuC7giHIuuhyAmHI+uhyAlJ1l3gUOSEHZqAVPWrY8KJyIuB1+xKjoSEhISEMCY5iSIid/X+/XoOx6yTkJCQcGIw1SLwT5srphX4EPD0yHgv2pVAO8ChyHoocsLhyHoockKSdRc4FDmnWQROSEhISJge6TKShISEhFOKpAASEhISTilmqQBE5LEi8tcicrWI/HDAX0Tk543/lSLyJTOV8ylGvitF5C0i8oAp5DSyDMrqhftSESlE5En7lM/Lf1ROEXmkiLxLRK4SkTftW0ZPjrH6v1BEXi0i7zayPnUiOV8qIteLSHCzxVz6k5FlTNZZ9KkxOb1wk/anUajqrP6AHPhb4HOBM8C7gS9ohXk88MdUH6h5KPBXM5Xz4cDtzPPjppAzVlYv3BuBy4EnzVFO4CLgfcA9zP93mmuZAj8C/JR5viPwceDMBLJ+BfAlwHt7/CfvTyvIOpc+NSin10Ym608xf3OcATwEuFpVr1HV81RfIru0FeZS4De0whXARa2tpbOQU1XfovUJ5yuAi/cso0VMmQI8C/h94Pp9CuchRs4nA69Q1b8DUNU5y6rAbUREgFtTKYDlfsUEVX2zybsPc+hPwLisc+lTEWUK0/enUcxRAdwN+Ij3/7XGbdUwu8aqMnwP1ShrCozKKiJ3ozqT8cI9ytVGTJl+HnA7EfkzEXmHiHzn3qRrIkbWXwQ+H7gOeA/wfaqz/AjAHPrTOpiyTw1iJv1pFHO8DC50S1N7r2pMmF0jWgYRuYSqsT5ipxL1I0bWFwA/pKqFTHdRVoycC+BBwFcBtwD+UkSuUNW/2bVwLcTI+hjgXcCjgHsDrxeRP1fVT+5YtlUxh/60EmbQp8bwAqbvT6OYowK4Fri79//FVCOoVcPsGlEyiMj9gZcAj1PVf9qTbG3EyPpg4OWmsd4BeLyILFX1VXuRsEJs3d+g1adKPy0ibwYeAOxbAcTI+lTgJ7UyCF8tIh8E7gu8dT8iRmMO/SkaM+lTY5hDfxrH1IsQgYWTBXANcC/qxbUvbIV5As1Fq7fOVM57AFcDD597mbbC/zrTLALHlOnnA/+PCXtLqmtE7jdTWX8ZeK55vjPwUeAOE7WBe9K/sDp5f1pB1ln0qTE5W+Em6U8xf7ObAajqUkSeCbyOahX9pap6lYg8w/i/kGpV/fFUDeEmqpHWHOX8MeCzgP9pRgJLneCWwEhZJ0eMnKr6fhF5LXAlUAIvUdW93yUVWabPA35dRN5DRa4/pKp7vyZYRP4P8EjgDiJyLfAc4MiTc/L+ZBEh6yz6VIScB4F0FURCQkLCKcUcdwElJCQkJOwBSQEkJCQknFIkBZCQkJBwSpEUQEJCQsIpRVIACQkJCacUSQEkJCQknFIkBZCQkJBwSvH/AyEiIJ4ZM5vuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make prediction\n",
    "prediction = model.predict(X_star, operator=None)\n",
    "\n",
    "u = griddata(X_star, prediction[:, 0], (X, T), method=\"cubic\")\n",
    "v = griddata(X_star, prediction[:, 1], (X, T), method=\"cubic\")\n",
    "\n",
    "h = np.sqrt(u ** 2 + v ** 2)\n",
    "\n",
    "\n",
    "# Plot predictions\n",
    "fig, ax = plt.subplots(3)\n",
    "\n",
    "ax[0].set_title(\"Results\")\n",
    "ax[0].set_ylabel(\"Real part\")\n",
    "ax[0].imshow(\n",
    "    u.T,\n",
    "    interpolation=\"nearest\",\n",
    "    cmap=\"viridis\",\n",
    "    extent=[t_lower, t_upper, x_lower, x_upper],\n",
    "    origin=\"lower\",\n",
    "    aspect=\"auto\",\n",
    ")\n",
    "ax[1].set_ylabel(\"Imaginary part\")\n",
    "ax[1].imshow(\n",
    "    v.T,\n",
    "    interpolation=\"nearest\",\n",
    "    cmap=\"viridis\",\n",
    "    extent=[t_lower, t_upper, x_lower, x_upper],\n",
    "    origin=\"lower\",\n",
    "    aspect=\"auto\",\n",
    ")\n",
    "ax[2].set_ylabel(\"Amplitude\")\n",
    "ax[2].imshow(\n",
    "    h.T,\n",
    "    interpolation=\"nearest\",\n",
    "    cmap=\"viridis\",\n",
    "    extent=[t_lower, t_upper, x_lower, x_upper],\n",
    "    origin=\"lower\",\n",
    "    aspect=\"auto\",\n",
    ")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Gxlw5jxzXD-L"
   },
   "source": [
    "**Assess the accuracy**\n",
    "\n",
    "We can employ the test data to assess the accuracy of the trained model.  \n",
    "Using the normalized $L^2$ norm defined as  \n",
    "  \n",
    "Error $u$ = $\\frac{1}{||u_{test}||_{L^2}} ||u_{test} - u_{pred}||_{L^2}$  \n",
    "  \n",
    "where $u_{test}$ is the reference solution and $u_{pred}$ is the prediction given by the model, we get:  \n",
    "  \n",
    "Error $u$: 1.854433e-03  \n",
    "Error $v$: 2.413796e-03  \n",
    "Error $h$: 1.426797e-03  \n",
    "  \n",
    "We can also plot the absolute value of the error for each point of the domain. It's everywhere in the order E-3.  \n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "j2SxMTchnP6C"
   },
   "source": [
    "![Error.png]()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "collapsed_sections": [],
   "name": "Schrodinger.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "py3.8",
   "language": "python",
   "name": "py3.8"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
